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Abstract 

This report of work performed under cooperative agreement NCC 2-755 with NASA Ames 
Research Center during the period August 1993 to January 1995 deals with the direct numerical 
simulation of transitional and turbulent flow at low Mach numbers using high-order-accurate 
finite-difference techniques. A computation of transition to turbulence of the spatially-evolving 
boundary layer on a heated flat plate in the presence of relatively high freestream turbulence was 
performed. The geometry and flow conditions were chosen to match earlier experiments. The 
development of the momentum and thermal boundary layers was documented. Velocity and tem- 
perature profiles, as well as distributions of skin friction, surface heat transfer rate, Reynolds shear 
stress, and turbulent heat flux were shown to compare well with experiment. The results indicate 
that the essential features of the transition process have been captured. The numerical method 
used here can be applied to complex geometries in a straightforward manner. 

Introduction 

Significant progress has been made in the past few years towards computing the flow fields 
associated with a variety of complex configurations. This can be attributed both to the increased 
speed and memory capabilities of computer hardware, and to improvements in computational 
algorithms and techniques. Complete three-dimensional solutions, both steady as well as 
unsteady, of the Reynolds-averaged Navier-Stokes (RANS) equations are possible on computers 
that are available currently. 

Despite these advances, the turbulence modeling that is required in RANS computations 
remains a major stumbling block. The quest for an accurate and universal turbulence model so far 
has proven unsuccessful. This is because turbulence models for the RANS equations are designed 
to model all scales of turbulent structures. Since the larger scales are geometry- and flow-depen- 
dent, it is possible that a universal turbulence model accounting for all scales of motion may never 
be developed. In addition to turbulence modeling, the inability to accurately predict flow transi- 
tion represents a second stumbling block. The solution to these problems might lie eventually in 
the use of techniques such as direct numerical simulations (DNS), where all pertinent scales of 
turbulence are computed using sufficiently fine grids that obviate the need for any modeling, and 
large-eddy simulations (LES), where only the smaller and more universal scales are modeled. 
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Despite their considerable computing requirements, a much greater reliance on such techniques is 
anticipated in the future as computer hardware capabilities improve. 

Most of the LES and DNS research to date has been performed using spectral-based methods. 
Further, research efforts have dealt primarily with the temporal evolution of transition in rela- 
tively simple geometries, and have been limited to low Reynolds numbers. Such temporal-evolu- 
tion” computations are approximations to the actual spatial evolution. In recent years, finite- 
difference DNS techniques are being explored as an alternative to spectral- based methods. Finite- 
difference methods have an advantage over spectral-based methods in terms of ease of applicabil- 
ity to complex geometries and implementation of general boundary conditions (required for spa- 
tially-evolving flows). Recent work 1,2 has shown that the accuracy levels required for DNS and 
LES can be achieved by using high-order upwind-biased finite-differences. Finite difference tech- 
niques appear well-suited for conducting DNS or LES of flows in complex and more realistic con- 
figurations. 

This report deals with a direct simulation of transition and turbulence in a spatially evolving 
boundary layer performed using high-order-accurate upwind-biased finite-difference techniques. 
The problem considered is the flow of a low Mach number fluid over a heated flat plate under rel- 
atively high freestream turbulence conditions. The geometry and flow conditions are chosen to 
match experiments 3 conducted at United Technologies Research Center and at NASA Lewis 
Research Center/Case Western Reserve University . 4 ’^’ 6 ’ 7,8 Our interest in the high freestream tur- 
bulence environment is two-fold. First, this is typical of several engineering applications, includ- 
ing gas turbine engines, that we hope to address in the future using the techniques developed here. 
Second, since the highly disturbed environment leads to an earlier onset of transition than in the 
low disturbance case, a much smaller computational domain (in the streamwise direction) is 
required. This in turn reduces drastically the computational effort. 

Boundary layer transition in a low-disturbance environment is characterized by a well-docu- 
mented sequence of events, starting with the initial amplification of two-dimensional Tollmein- 
Schlichting waves, followed by the formation of three-dimensional spanwise-periodic structures, 
development of turbulent spots, and finally culminating in random turbulent motion. The methods 
of linear stability theory have been used successfully to predict the Reynolds number above which 
small linear disturbances are amplified. However, in highly disturbed environments the presence 
of large disturbances may lead directly to turbulent spot formation. Morkovin^ has used the term 
“bypass transition” to describe this situation where the linear instability mechanisms are bypassed 
by finite-amplitude nonlinear instabilities. This latter mode of transition is not very well under- 
stood. 

One of the main motivations for the present work is to improve our fundamental understand- 
ing of thermal and momentum transport in transitional flows in the presence of high freestream 
turbulence. Such an understanding is important for many engineering applications. It is of 
extreme importance in improving heat transfer modeling and thermal load predictions, and deter- 
mining cooling requirements for aircraft engine gas turbine blades-especially in the first stage of 
the high pressure turbine which is influenced largely by the combustor exit turbulence and where 
the flow over a large portion of the stator airfoil is typically in a transitional state. Current design 
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practice for turbine blading compensates for our lack of understanding by overcooling, resulting 
in aerodynamic cooling penalties, lowered burner and turbine mainstream mass flow rates 
(reduced power/efficiency), and potentially increased cooling system fabrication costs . The 
present work eventually could impact turbine airfoil design by providing both a fundamental 
understanding of the flow phenomena, and information for improved modeling/computation of 
boundary layer transition in turbomachinery. To this end, central to this work is the creation of an 
extensive data base that can be used for future modeling studies, as well as validation of RANS 
and LES computations. 

An important aspect of the present work is that we hope to shed some light on certain unex- 
pected results that have been encountered by researchers in measuring the turbulent normal heat 
flux, V? , in the transitional boundary layer on a heated flat plate using three- wire probes. In the 
transition region, negative values were reported 0 near the wall for this correlation between the 
fluctuating normal velocity and the fluctuating temperature. This is an unexpected and counter- 
intuitive result since it implies a counter-gradient heat transfer by the fluctuations. Since the mean 
temperature gradient is negative, a negative value for the time-averaged v t implies a negative 
eddy thermal diffusivity which would appear improbable. More recent measurements 8 in the 
same facility but using a different set of probes have duplicated this result, as have measurements 
reported by other researchers in a completely different facility. 10,11,12 Various reasons for the 
negative v't' correlation, including spatial resolution of the probe, wind-tunnel peculiarities, 
insufficient frequency response of the cold wire, thermal crosstalk between the hot and cold wires 
and near-wall “streaky structures”, have been suggested and some of these have been subse- 
quently disproved by measurements. 8,10 In contrast, measurements at the University of Minne- 
sota for both flat and concave walls ^ do not exhibit negative values for the v't' correlation. 

An additional aspect of the present work is the evaluation of the turbulent Prandtl number in 
the transition region to aid in modeling heat transfer from turbine blades. Measurements for this 
quantity in turbulent boundary layers have been reported by several researchers 1 ^’ 1 ^’ 10 but such 
data is lacking for transitional boundary layers. 8 This is because boundary layers in highly dis- 
turbed freestream environments undergo transition close to the leading edge, resulting in thin 
boundary layers that make accurate measurement somewhat difficult. The only direct measure- 
ments of the turbulent Prandtl number in transitional flows were made by Kim et al. 17 who 
reported measurement uncertainties higher than 20%. Other experiment al ef forts aimed at mea- 
suring this quantity 7,11 were beset by problems relating to the negative v't' measurement men- 
tioned above. 

The present work complements earlier DNS computations of boundary layer transition on a 
flat plate under adiabatic conditions. 2 A direct comparison can thus be made between the present 
results for a heated wall and these earlier results for an adiabatic wall. One aspect in which the 
present work differs from this earlier work 2 is in the treatment of the leading edge of the flat plate. 
Earlier, a sharp leading edge was considered; the degree to which the leading edge singularity 
affected the downstream flow remained a matter of conjecture. Here, the elliptical leading edge of 
the plate is included in the computations. This will enable leading edge receptivity studies to be 
performed. 
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Geometry and Grid System 

The computational domain that was used was modeled after the configurations of prior exper- 
iments conducted at United Technologies Research Center^ and at NASA Lewis Research Center/ 
Case Western Reserve University . 4 ’ 5,6,7 ’ 8 Both sets of experiments are essentially similar, differ- 
ing slightly in the length of the unheated starting length and in the input heat flux. The unit Rey- 
nolds number based on inlet conditions is 50,000/inch, and the inlet Mach number, M , is 0.09. 
The freestream turbulence intensity in the experiments is approximately 2.6%, the onset of transi- 
tion is at roughly Re x = 250,000, and the flow becomes turbulent around R& x = 600,000. The flat 
plate used in the computation has a 4: 1 elliptical leading edge. As in the experiment^, an unheated 
starting length of 1.7 in. (including the elliptical leading edge) is provided upstream, beyond 
which the flat plate is maintained at a constant heat flux of 0.075 Btu/sec.ft 2 . 

The computational region of interest (Fig. 1) consists of several zones. Zone 1 is an inlet 
region where numerically generated freestream turbulence develops. Zone 1 exchanges informa- 
tion with three different zones: zone 2, which covers the lower portion of the plate, zone 3, which 
covers the elliptical leading edge of the plate, and zone 4 which spans the region above the plate. 
Zone 4 represents the primary region of interest here and contains the highest grid resolution. 
Downstream of zone 4 is an exit region (zone 5) where the grid becomes gradually very coarse in 
the streamwise direction. 


10 in. ^ 

15 in. 

12 in. 


Inlet Zone 

Well Resolved Region of Flow 

Exit Zone 

a) 

( 4 ) 

( 5 ) 

Leading Edge 
Zone (3) 

Plate Lower Zone (2) 


Figure 1 . Flow geometry and zonal grid 
configuration used in the direct simulation. (Not 
to scale.) 

The streamwise extent of the inlet region (zone 1) is 10 inches. The well-resolved region of 
interest on the flat plate (zone 4) spans a total of 15 in., and the exit region (zone 5) extends 
another 12 inches. Thus, the well-resolved region in the computation extends downstream of the 
plate to Re v = 750,000, and includes the transition region and a portion of the turbulent region on 

* X ... 1 8 

the plate. The momentum thickness at this location is 0.031 in. (based on the expression 
Q/x = 0.0142R<?~ I/7 ), or Re Q = 1540. The spanwise dimension of the computational region 
was chosen as 1.6 in., and is similar to that used in previous work. 2 The transverse extent (height) 
of the computational domain is 3.0 in. above and below the flat plate. (The height of the inlet 
region is thus 6.0 in.). Since symmetry boundary conditions are used on the upper and lower 
boundaries of all the zones, the simulation is in effect a channel-flow simulation, albeit a very 
wide channel (approximately 100 times the momentum thickness at Re x = 750,000). The result- 
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ing acceleration of the flow is extremely small and can be neglected for purposes of comparison 
with experiment. 

The grid system consisted of a 102 X 55 x 61 grid in zone 1, a 35 x 55 X 61 grid in zone 2, a 55 
x 40 x 361 curvilinear grid in zone 3, a 1030 x 55 x 361 in zone 4, and a 97 x 55 x 361 grid in the 
exit zone 5. Stretched grids are used in the transverse direction with a minimum Ay = 1 , based 
on the wall shear velocity at Re x = 750,000. Uniform grids are used in the spanwise direction, 
with Az + = 10. For the most part, uniform grids are used in the streamwise direction in the main 
zone (zone 4), with Ax + = 27. Stretching in the streamwise direction is used in the initial portion 
of zone 4 where the flow is laminar, in zone 5 which forms the exit region, and in zone 2. The total 
number of grid points used is roughly 23.63 million (of which 22.38 million are in zone 4). 

Numerical Method 

In the present approach, the time-dependent compressible Navier-Stokes equations are solved 
in nonconservative form in three-dimensional Cartesian coordinates (except in the region around 
the leading edge where curvilinear coordinates are used). High-order accurate, upwind-biased 
finite-differences are used for spatial discretization and the solution is advanced in time using an 
iterative-implicit scheme. Upwind-biased differences are used here because they achieve higher 
orders of accuracy on compact stencils than fully upwind schemes. The loss of accuracy near 
computational boundaries where the stencil size available is smaller is thus limited. The use of the 
nonconservative form of the Navier-Stokes equations limits the present method to flows that are 
free of discontinuities (such as purely subsonic flows). 

The convective terms are computed using fifth-order-accurate forward- and backward-biased 
finite-differences on a seven-point stencil. Viscous terms are computed using fourth-order-accu- 
rate central differences and a five-point stencil. Since we are dealing with stretched meshes, the 
coefficients in the difference formulas are evaluated numerically using Lagrange polynomials so 
as to retain high-order accuracy even on grids where the rate of change of grid spacing is not suf- 
ficiently smooth. Obtaining explicit formulas for these coefficients is a laborious task for differ- 
ences with higher than second-order accuracy. Additional details regarding this procedure, and the 
evaluation of the convective terms and viscous terms, can be found in previous work 2 . 

The solution is advanced in time using an iterative-implicit approach where the fully implicit 
form of the equations (not the linearized implicit equations) are solved at each time step. Factor- 
ization and linearization errors are thus driven to zero at each time step. The overall method is 
thus fourth-order accurate in space, second-order accurate in time, and does not require the speci- 
fication of any arbitrary smoothing parameters. 

Boundary Conditions 

The conditions imposed at the various boundaries of the computational domain are discussed 
here. The inlet boundary of the computational domain (left boundary of zone 1 in Fig. 1) is a sub- 
sonic inlet boundary where four quantities must be specified. The four chosen for this study are a 
Riemann invariant R j = u + 2c/ (y- 1) , the stagnation pressure, and the velocities in the y and 
z directions, v, w. Here u denotes the streamwise velocity and c the local speed of sound. A sec- 
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ond Riemann invariant, R 2 = u-2c/(y- 1) is extrapolated from the interior of the domain to 
complete the system of equations required to uniquely determine the dependent variables at the 
inlet. As was done in Rai and Moin, 2 small velocity perturbations are introduced at the inlet 
boundary. These perturbations then evolve through the domain to yield a disturbance field at the 
exit to the inlet region (i.e., leading edge of the plate) that closely resembles the experimental 
freestream turbulence. Quantities that are matched with the experiment include the turbulence 
intensities, the longitudinal integral length scale, and the power spectrum. Further details can be 
found in earlier work. 2 

The exit boundary of the domain (right boundary of zones 5 and 2) is a subsonic exit boundary 
where one quantity must be specified. Here the exit static pressure is specified and the remaining 
variables are extrapolated from the interior of the domain. The advantage of using this boundary 
condition in conjunction with the inlet stagnation pressure specification is that the mass flow 
through the domain is uniquely determined. However, this is a reflective boundary condition that 
reflects pressure waves back into the domain. The effects of these reflections are mitigated by 
gradually coarsening the grid in the exit region of the domain in the streamwise direction (by a 
factor of about 100). 

The coarsening of the grid in the exit region serves an additional purpose. The flow across the 
exit boundary of the computational domain consists of an inviscid outer region (which comprises 
most of the area of this boundary) and the boundary layer. The static pressure from the inviscid 
region (in the vicinity of the boundary layer) is imposed within the boundary-layer region of the 
exit boundary. However, this approach is valid only when the pressure gradient normal to the wall 
is small (and is not the case in an unsteady turbulent boundary layer). The grid coarsening results 
in the unsteadiness within the boundary layer being dissipated numerically. 

Symmetry boundary conditions are imposed on the upper and lower boundaries of the 
domain. This is done by creating temporary arrays of variables above (or below) the plane of sym- 
metry such that the equations of motion can be integrated at grid points on the symmetry plane. 
Dependent variables above the plane are computed by reflection. Periodicity conditions are 
imposed at the spanwise boundaries of the domain. 

Along the flat plate, the “no-slip” boundary condition is used together with a constant heat 
flux condition (or adiabatic wall in the unheated starting length region) and zero normal pressure 
derivative condition. 

The zonal boundaries that separate the various zones are treated by using a fourth-order accu- 
rate interpolation procedure that is explicit at each iteration within a given time step, but implicit 
over the entire time step 2 (The transfer of information between the inlet zone (zone 1) and the 
leading edge zone (zone 3) is, however, only second-order accurate.) Owing to computer memory 
constraints, the main zone of interest (zone 4) is divided up into several blocks and the governing 
equations are solved in each block in a sequential manner. This introduces additional artificial 
zonal boundaries. Data along several streamwise planes at either ends of each block are stored in 
additional arrays that are accessed when evaluating the fluxes thus maintaining the same spatial 
accuracy at these boundaries as at interior points. The only approximation introduced by these 
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artificial zonal boundaries is in the linearization, but this is driven to zero by the multiple itera- 
tions that are performed at each time step. 


Results 

The numerical method and boundary conditions discussed in the previous sections and the 
corresponding computer program were first validated by studying the growth of two-dimensional, 
small amplitude disturbances in flat plate flow. This investigation closely resembled that of Fasel 
et al. 19 A similar investigation was carried out to validate the computer code in Rai and Moin. In 
the present investigation the inlet region, the plate leading edge, and the zonal boundaries are 
included and their effect on the results were documented. The results of these validation studies 
are not presented here, but they show excellent agreement with earlier work * and with linear 
theory. 

Results obtained from the direct simulation of boundary layer transition on a heated flat plate 
are presented in this section. These results were computed by integrating the governing equations 
in conjunction with the boundary conditions (both natural and zonal) described in earlier sections. 
Three iterations were performed at each time step. The root-mean-square residual decreased by a 
factor of about 15 at the end of the third iteration in the main zone of interest (zone 4). Although a 
larger decrease in the residual is desirable, the number of iterations was restricted to three both in 
the interest of obtaining a solution in a reasonable amount of time and because detailed studies 
reported earlier 2 noted that three iterations provided adequate solution accuracy. 

The results presented here were obtained in a series of stages. The initial solution was 
obtained on an extremely coarse grid and then interpolated onto successively finer grids. The 
computation on each grid was stopped when the skin-friction distribution (after spanwise averag- 
ing) showed minor changes in time and exhibited an average in time that essentially was in equi- 
librium. 

The nondimensional time-step used in the computation was Atu^/8 = 0.04, where 5 is 
the boundary layer displacement thickness at the end of the well-resolved region of the plate, Re x 
= 750,000. The maximum CFL number corresponding to this time-step is about 40.0 (the CFL 
number based on the mean convective velocity is approximately 4.0). 

The results presented here are based on a rather small statistical sample roughly 5 millisec- 
onds in duration. This sample size is marginally adequate for the results shown here. Much of the 
jitter in the results presented here is a reflection of the small sample size. More detailed analyses 
of the flow field such as conditional sampling or intermittency studies would require a sample that 
is at least ten times larger. The computer program requires about 4.5 microseconds of single-pro- 
cessor Cray C90 time per grid point per iteration. Approximately 300 hours of Cray C90 single- 
processor time was required to acquire the current sample. 

Characterization of Freestream Disturbance 

As described in the Boundary Conditions section, the freestream disturbance is generated by 
perturbing the three velocity components at the inlet of zone 1 to result in freestream disturbances 
at the exit to zone 1 (which is the plate leading edge) that matched the experimental conditions at 
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this location. The inlet perturbation velocities are obtained using a prescribed power spectrum. 
The longitudinal integral length scale, was taken to be 0.5 in. The intensity values at the inlet 
of zone 1 were iteratively modified until the intensity values at the leading edge of the plate were 
approximately the same as those obtained in the experiments^ (where it is referred to as the grid 
No. 2 case). 

Figure 2 shows the variation of u\ v', w', and the freestream turbulence intensity T rm in the 
inlet region as a function of the streamwise distance. Within a few inches from the inlet to zone 1, 
a decaying, nearly isotropic turbulence field is set up and the freestream turbulence intensity T rms 
is seen to match well with the theoretical -5/7 decay law (curve marked theory in Fig. 2) for grid 
generated turbulence. At the exit of this region (which is the leading edge of the plate), the turbu- 
lence intensity is about 2.6%, which is the value reported in the experiments. 3 There is some vari- 
ation in the intensity values for the three components. Note that the experimental data 3 also show 
a variation of ±15% in intensity values at the leading edge with v' > w' > u . 



Figure 2. Variation of free-stream turbulence 
intensities in zone 1 . 


Figure 3 shows the power spectrum obtained immediately upstream of the leading edge of the 
flat plate for the streamwise component of the velocity. This power spectrum was computed using 
fast Fourier transforms in conjunction with a Hanning window. 20 The size of the window was 
chosen as the period of the lowest frequency perturbation introduced at the inlet to zone 1 (the 
lowest frequency being 24Hz). The computed spectrum roughly approximates the Von Karman 
spectrum. We have performed more detailed investigations that have yielded a better understand- 
ing of what is required to accurately reproduce experimentally observed freestream turbulence. 
This work was done under the present cooperative agreement during the period March 1992- July 
1993 and is not discussed here. 


9 





Figure 3. Computed power spectrum for the streamwise 
velocity component at streamwise location immediately 
upstream of the leading edge of the plate. 


Surface Skin Friction and Heat Transfer 

Figure 4 shows the computed skin friction along with experimental data, 4 ’^ the Blasius solu- 
tion for laminar flow, and the turbulent correlation^ for the skin friction distribution. The Blasius 
solution is given as: 

C, = 0.664 Re~ l/l ( 1 ) 

The turbulent correlation is given as : 

C f = 0.455//n 2 (0.06/?^) (2) 

The computed data indicates onset of transition at roughly the same location as the experi- 
ment. In the transition region the computed skin friction increases rapidly and then follows the 
turbulent correlation. The computed rise in skin friction is much more rapid than the experimental 
measurement. However, there is considerable uncertainty associated with the experimental skin 
friction values since these are derived quantities obtained from the mean velocity profiles using 
the momentum integral relation. The presence of small oscillations in the computed skin friction 
is a reflection of the small sample size. Smoother data can be obtained by using a larger sample 
size. 
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■ Sohn et al. 
o Present Results 



Figure 4. Computed skin-friction coefficient 
distribution along the length of the plate. 

o £ 

Figure 5 shows the computed Stanton number variation. Two sets of experimental data ’ are 
also shown in the figure, along with the expressions for laminar and turbulent boundary layers 
with corrections to account for the unheated starting length/*’ 2 * For laminar boundary layers, the 
Stanton number is given as: 

St = 0.453Pr -2/3 fl<?; 1/2 (l - (x 0 /x) 3/4 ) ^ (3) 

where x Q represents the unheated starting length and Pr is the Prandtl number (= 0.72). For 
turbulent boundary layers, the Stanton number is given as: 

St = 0.0 3Pr~ 0A Re~ 02 ( 1 - (x/x) 0 ' 9 )^^ (4) 

The computed data agree well with the laminar variation. The computed onset of transition 
occurs slightly further downstream than in the experiments. Overall, there is reasonable agree- 
ment between the computation and experiments, with the computed results being slightly lower in 
the transition and turbulent flow regions. The agreement is clearly much better than for the skin 
friction (Fig. 4) and can be attributed largely to the higher accuracy of the Stanton number mea- 
surement. The two sets of experimental data shown in the Figure are in close agreement with each 
other in the transition and turbulent regions. In the laminar region, the differences between the 
two experimental data sets is because of slightly different unheated starting lengths. The effect of 
unheated starting length is much greater for laminar than for turbulent boundary layers. While the 
computed skin friction data seemed to agree qualitatively with the turbulent correlation beyond 
transition, the level of the Stanton number curve is slightly lower than the turbulent correlation. 
This suggests that higher grid resolution is required to match the Stanton number than required to 
match the skin friction. 
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Figure 5. Computed Stanton number distribution 
along the length of the plate. 


Reynolds Analogy Factor 

A plot of the Reynolds analogy factor, 2 St/ Cp is shown in Fig. 6. The Reynolds analogy fac- 
tor in air for a flat plate with zero pressure gradient and a constant wall temperature boundary con- 
dition is well represented by Pr~ 2/21 for both laminar and turbulent flows. However, the reference 
curves shown in the Figure for the laminar and turbulent regimes differ from the conventional 
Pr~ 2/2> curve due to the effects of the uniform wall heat flux boundary condition and the 
unheated starting length. The curves shown were obtained by combining appropriate laminar and 
turbulent theoretical results. 6,21 The laminar variation plotted in the Figure is given as: 

2St/C f = 1.364/V _2/3 (1 - (x 0 /x) 3/4 ) W * (5) 

The turbulent variation is expressed as: 

—0 4 09 -1/9 

2 St/C f = 1.045 Pr U ' ( 1 - ( *,/*) ) (6) 

The computed data in the laminar and turbulent region agree quite well with the correspond- 
ing theoretical curves. The computed data in the early transition region agree with the experiment, 
but dip below the turbulent correlation in the late transition region before approaching the fully 
turbulent value. This behavior is a result of the rapid increase in skin friction that was noted in 
Fig. 4. 
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Figure 6. Computed Reynolds analogy factor 
distribution along the length of the plate. 


Boundary Layer Thickness 

The development of the boundary layer thickness in the streamwise direction is shown in 
Fig. 7. The local boundary layer thicknesses shown here were used to nondimensionalize the dis- 
tance in the wall-normal direction in several of the figures shown later in this report. Also 
included in Fig. 7 are the laminar boundary layer thickness from the Blasius solution and the tur- 
bulent boundary layer thickness obtained from a momentum-integral analysis for a turbulent 
1 ft 

boundary layer. 
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thickness along the length of the plate. 

The laminar boundary layer thickness is obtained from the relation: 

Re 5 = 5X)Re\ /2 

The turbulent boundary layer thickness is obtained from the relation: 


(7) 


Re 5 = 0.14 Re 6 / 1 (8) 

The computed boundary layer thickness agrees with the laminar variation and gradually 
increases in the transition region to match the turbulent values. 

Mean Velocity and Temperature Profiles 

Figure 8 shows computed mean velocity profiles obtained at various streamwise locations. 
The dashed lines represent the linear relationship u + = y + and the logarithmic law of the wall 
u + = 2.5 lny + + 5.0. The laminar or Blasius velocity profile at Re x = 250,000 is also shown. The 


14 




computed velocities gradually change from laminar to the turbulent profile with increasing Re . 



Figure 8. Mean velocity profiles at various 
streamwise locations, in wall units. 


The changes in the velocity profiles are more rapid in the transition region than in the early turbu- 
lent region. Similar trends can be seen in the experimental data. 3,8 In the experiments, turbulent 
profiles are attained earlier in one experiment 8 than in the other. 3 The computed profiles indicate a 
transition region that lies between the two sets of experimental data. In Fig. 8, the profiles in the 
laminar region close to the plate leading edge show a dip in the outer portion of the boundary 
layer. The reasons for this are currently being investigated and appear to be related to leading edge 
effects. The dip gradually disappears as the streamwise distance from the leading edge of the plate 
increases. The computed turbulent profile at Re x = 725,000 satisfies the law of the wall and the 
log-law. A significant wake region is not seen in the profile. Although the experiments clearly 
showed that the wake strength was diminished with increasing freestream turbulence, a small 
wake region was noted in the experiments at the present freestream turbulence level. The reason 
for this difference is not clear. 

Figure 9 shows computed mean temperature profiles obtained at various streamwise locations. 
The dashed lines represent the equation t + = Pr. y + and the temperature log-law curve for fully 
turbulent boundary layers 6 : 


' + = l3 ' 2Pr+ m ln &2 


) 


(9) 


where Pr is the Prandtl number (= 0.72) and Pr ( is the turbulent Prandtl number, assumed to 
be a constant value of 0.9. The temperature in wall units, t + , is defined as t + = (T —T) /f x , 
where t x is referred to as the friction temperature and is related to the wall heat flux, q w , and the 
friction velocity, « x , as t x = q w /pC u v The agreement between the temperature profile in the 
turbulent region Re x = 725,000 and the temperature log-law is approximately the same as shown 
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in Fig. 8 for the mean velocity profiles, although a wake region can be discerned in Fig. 9. A com- 
parison of the profiles in Fig. 9 and Fig. 8 do not show any marked lag (or lead) in the develop- 
ment of the mean temperature compared to the mean velocity profiles. Several experimental 
studies 10,17,22 have noted a lag in the development of the mean temperature compared to the 
mean velocity profiles, but the experiments of Sohn et al. 7 show an opposite trend. 



Figure 9. Mean temperature profiles at various 
streamwise locations, in wall units. 


RMS Velocity and Temperature Profiles 

Figures 10 and 1 1 show the evolution of the streamwise component u' of the turbulence inten- 
sity normalized by the freestream velocity at various locations along the plate. In Fig. 10 wall 
units are used, while in Fig. 1 1 outer variables are used. Both figures show that the peak value of 
u' increases rapidly with increasing Re x and reaches its maximum (about 18.5% of the free- 
stream velocity) at about Re x = 400,000 in the transitional region. Thereafter it decreases to about 
15% of the free-stream velocity in the turbulent region. In terms of wall units, in the laminar and 
pretransitional region of the flow (Re = 200,000 and 300,000), the peak value is located at 
approximately y =27; the peak shifts to a location y = 17 in the transition region and 
remains at this location through the early turbulent region (Fig. 10). In terms of actual distance 
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from the wall, the location of the peak value can be seen moving closer to the wall with increasing 
Re x in Fig. 11. 



Figure 10. Streamwise component of turbulence 
intensity (normalized by ffeestream velocity) at 
various streamwise locations, plotted in wall 
coordinates. 
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Figure 11. Streamwise component of turbulence 
intensity (normalized by freestream velocity) at 
various streamwise locations. 

The general trend of the profiles in both Fig. 10 and Fig. 1 1 is the same as that observed exper- 
imentally. 6,4,10 However, one feature that is seen in the experiments and not in the computation is 
the presence of a second peak in each of the profiles in the transition region. The second peak is an 
artifact of the experimental technique and occurs because of the switching between laminar and 
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turbulent boundary layer flows as a turbulent spot goes by the hot-wire probe. The switching gives 
rise to a velocity fluctuation level that is higher than that obtained in either the laminar or the tur- 
bulent regions and affects the profile. 
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Figure 1 2. Normal component of turbulence 
intensity (normalized by freestream velocity) at 
various streamwise locations. 



Figure 13. Spanwise component of turbulence 
intensity (normalized by freestream velocity) at 
various streamwise locations. 


Figure 12 shows the normal component v' of the turbulence intensity at the same axial loca- 
tions as the u' profiles shown in Fig. 11. A considerably different evolution is observed in the v' 
profiles than in the u' profiles. The peak value increases rapidly in the transition region with 
increasing Re to a maximum of about 5% and then maintains this value through the early turbu- 
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lent region. No significant decrease is observed to occur in v' once the maximum value is 
attained, unlike the behavior of u' which decays in the late transition/early turbulent stage. Pro- 
files of the spanwise component w' of the turbulence intensity are shown in Fig. 13. The general 
trend in these profiles is similar to that seen in the v' profiles. Figure 14 shows the variation of the 
peak value of the three turbulent intensity components with Re . The trends in the peak values 
discussed earlier are clearly seen in this figure. 



Figure 1 4. Variation of peak turbulence intensities 
with streamwise distance. 


The computed turbulence intensities in the early turbulent region of the flow (at Re x = 
725,000) are compared to experimental data in Fig. 15. The symbols used in this Figure represent 
experimental data 24 obtained in a closed water channel using the LDV technique. The experimen- 
tal u' and v' data were obtained at a Reynolds number based on momentum thickness, Re§, of 
2420. The experimental w' data were obtained at Re^ = 1750. The computed data are at a stream- 
wise location corresponding to Re§ = 1350. The free-stream turbulence level in the experiment 
was approximately 1.8%. 

Except for u' in the region 8 < y + < 100, the computed intensities are in good agreement with 
the experimental data. The nature of the discrepancy (higher u' values and slightly lower v' val- 
ues) has been observed in earlier channel-flow simulations with coarse grids. 1 Additionally, com- 
putations on coarser grids in the present study clearly indicate that the major contributor to the 
discrepancy is the coarse grid used in the computation. The differences in the computed and 
experimental data for values of y + > 500 are mainly caused by the free-stream turbulence level in 
the experiment being different from that in the computation. Based on the coarse grid studies per- 
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formed here and on grid refinement studies performed earlier, 2 we expect the predicted peak value 
to drop to the experimental value at the next grid refinement level. 



Figure 15. Turbulence intensities at streamwise 
location Re x = 725,000; normalized by wall- 
shear velocity and plotted in wall coordinates. 



Figure 16. Turbulence intensities at streamwise 
location Re x = 725,000; normalized by local 
mean velocity and plotted in wall coordinates. 


The computed turbulence intensities normalized by the local mean streamwise velocity at the 
same streamwise location as in Fig. 15 ( Re x = 725,000) are shown in Fig. 16. The symbols repre- 
sent the experimental data. 24 The streamwise and spanwise components agree well with the 
experimental data. The limiting value at the wall of the streamwise component is very nearly the 
same as the experimental value of 0.39. The normal component is lower than that of the experi- 
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mental data. It should be noted that the experimental v'/u data do not exhibit the expected linear 
behavior in the near-wall region; instead, they increase as the wall is approached thus indicating 
experimental errors. 

The velocity profiles and intensity profiles presented are noted to be in good agreement with 
earlier computations 2 performed for an adiabatic flat plate under nearly identical flow conditions. 
This is to be expected, since the wall heat flux in the present study is small. 



Figure 17. rms temperature fluctuations at various 
streamwise locations, plotted in wall coordinates. 
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Figure 18. rms temperature fluctuations at various 
streamwise locations. 

The rms temperature fluctuations at various streamwise locations are shown in wall units in 
Fig. 17 and in outer variables in Fig. 18. In both figures the temperature fluctuations are scaled by 
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(7 - 7 ^) , where T w denotes the local mean wall temperature. In terms of wall units, the tem- 
perature fluctuations show a plateau region up to about y + values of 10 - 20 and then drop off in 
the outer portion of the boundary layer. In the early turbulent region, the drop is more gradual than 
in the transitional and laminar regions. The same data when plotted in outer variables (Fig. 18) 
show the migration of the peak toward the wall. Both figures show the intensity of the temperature 
fluctuations increasing rapidly through the transition region and then decreasing in the early tur- 
bulent region. 

Reynolds Shear Stress and Heat Flux Profiles 

Figures 19 and 20 show Reynolds shear stress profiles at the same streamwise locations as the 
u' profiles shown in Fig. 1 1 . Clearly, the peak values of these profiles follow the same trend as the 
u profiles in Fig. 1 1 . As expected, the peak values close to the leading edge are much smaller 
than those in the turbulent region. There is a rapid increase in peak value in the region 
200, 000 < Re x < 400, 000. Thereafter the peak value gradually decreases to the turbulent value 
at Re x = 725,000. Figure 20 shows a region extending from about y + = 30 to 200 in the early tur- 
bulent region where the turbulent shear stress is approximately constant. 

Figure 21 compares the computed Reynolds shear stress values in the early turbulent region, 
Re r = 725,000, with experimental data. 24 Overall there is good agreement between experiment 
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and computation. Also, the results are very similar to earlier computations for the adiabatic plate. 
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Figure 19. Reynolds shear- stress distributions at 
various streamwise locations; normalized by the 
square of the wall-shear velocity and plotted in 
outer variables. 
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Figure 20. Reynolds shear stress distributions at 
various streamwise locations; normalized by the 
square of the wall-shear velocity and plotted in 
wall coordinates. 



Figure 21. Reynolds shear stress distribution at 
the streamwise location Re x = 725,000; 
normalized by the square of the wall-shear 
velocity and plotted in wall coordinates. 


We now turn our attention to the turbulent heat flux distributions. As noted in the Introduction, 
one of the motivations for the present study was to shed some light on experimentally observed 
anomalies in the measurement of the turbulent heat flux. In particular, certain experiments showed 
unrealistic negative values for the turbulent heat flux near the wall in the transition region. These 
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unexpected experimental results still remain unresolved. Our purpose here is to determine the 
behavior of the turbulent heat flux through the transition region to determine if there is any basis 
to the experimental observations. Figures 22 and 23 show the turbulent heat flux distributions at 
various streamwise locations, plotted in outer and inner variables, respectively. In both the figures, 
the turbulent heat flux, vY, is normalized by q w /pC p . There is clearly no evidence of negative 
values for vY either in the transitional or in the turbulent region. Both Fig. 22 and Fig. 23 bear 
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Figure 22. Turbulent heat flux distributions at 
various streamwise locations; normalized by 
q w / p C n , and plotted in outer variables 
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Figure 23. Turbulent heat flux distributions at 
various streamwise locations; normalized by 
q w /pC n , and plotted in wall units. 
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close resemblance to the Reynolds shear stress distributions shown in Fig. 19 and Fig. 20, respec- 
tively. 

Figure 24 shows the turbulent heat flux data normalized by the rms values of v' and t', to 
allow comparisons with the limited data reported in the literature for turbulent boundary layers. 
Figure 24 shows that in the turbulent region of the flow ( Re x = 725,000), the turbulent heat flux 
values range from 0.4-0. 5 in the boundary layer. These values are in agreement with measure- 
ments reported in the literature for turbulent boundary layers, where values ranging from 0.42 to 
0.54 have been observed. 8 ’ 14, 6 Experimental data in the near-wall region is, however, lacking. 



Figure 24. Turbulent heat flux distributions at 
various streamwise locations; normalized by the 
rms values of v' and t'. 


Turbulent Prandtl Number 

Figure 25 shows the computed turbulent Prandtl number distribution at various streamwise 
locations along the plate. At all streamwise locations, the curves are seen to rise sharply away 
from the wall and then drop in the outer region of the boundary layer. The location of the peak 
value moves closer to the wall with increasing streamwise distance. As mentioned earlier in the 
Introduction, to the best of our knowledge, there is only one reported set of experimental data 17 
(with rather large uncertainties) for the direct measurement of Pr ( in transitional flows. These 
measurements also showed turbulent Prandtl number values (sampled on the intermittency func- 
tion) substantially below unity in the near- wall region, implying that the eddy diffusivity of heat 
increases during transition more rapidly than does the eddy diffusivity of momentum. The 
available experimental data for Pr , variation in turbulent boundary layers show striking inconsis- 
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tencies; measured data near the wall is also somewhat lacking. Some experimental results show 
Pr ( values of unity and slightly higher everywhere in the boundary layer) while other experi- 
ments 16 show Pr ( decreasing from a peak of unity to a value of about 0.5 at the edge of the 
boundary layer. Outside the near-wall region, the computed results in the turbulent portion ( Re x = 
725,000) seem to agree with the latter set of experiments. Note also that this behavior is similar to 
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the distribution suggested by Rotta. 25 . Near the wall, Pr ( drops sharply in both the transitional 
and turbulent region. This is in contrast to some experimental measurements 21 in turbulent 
boundary layers that show Pr rising sharply above unity near the wall. 



Figure 25. Computed turbulent Prandd number 
distribution at various streamwise locations along 
the plate. 


Instantaneous Flowfield Visualizations 

Some limited visualizations of the instantaneous flowfield were performed using vorticity 
contours. More detailed visualization efforts to extract various features of the flow are currently 
underway. 

Spanwise vorticity contours in an (x, z) plane at > ,+ =1.0 (based on the wall shear velocity at 
Re = 725,000) are shown in Fig. 26 and in Fig. 27. Figure 26 shows contours in the region 
300, 000 < Re x < 500, 000. Isolated patches of vorticity are seen immersed in relatively quiescent 
regions. Vortical structures are seen to appear rather abruptly in the vicinity of Re x = 375,000. 
Figure 27 shows contours in the region 500, 000 <fl<? x <700, 000. The transition process is 
essentially completed in this region. These plots of spanwise vorticity are quite similar to earlier 
results 2 presented for a flat plate under adiabatic wall conditions. This is to be expected, since the 
freestream conditions in both cases are almost identical, and the wall heating used here is quite 
small. 
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Figure 26. Spanwisc vorticily contours in an (v, ;) plane at y - 1.0 in the region 

300, 000 < Re x < 500. 000. 



Figure 27. Spanwisc vorticily contours in an (.x; ;) plane at y — 1.0 in the region 

500, 000 < Rc x < 700, 000. 




Summary 

A direct simulation of transition and turbulence in a spatially evolving boundary layer on a 
heated flat plate is described. The simulation is performed using a high-order-accurate, upwind- 
biased, iterative-implicit, finite-difference algorithm. The algorithm is developed for the unsteady, 
compressible form of the Navier-Stokes equations. However, the simulation is performed at a low 
subsonic Mach number (0.09) because experimental data are available for this speed regime. This 
study focuses on the high-freestream disturbance case in which transition to turbulence occurs 
close to the leading edge, thereby substantially reducing computing requirements. A zonal 
approach that allows efficient and judicious distribution of the grid points is used. The zonal pro- 
cedure also allowed the leading edge of the plate to be included easily in the computation. No vis- 
ible degradation of solution quality was observed in the vicinity of the zonal boundaries. 

A variety of information was extracted from the simulation in order to elucidate the nature of 
the computed transition process and provide comparisons with experimental observations. Some 
characteristics of the momentum and thermal boundary layers have been documented. The results 
indicate that the essential observed features of the transition process have been captured in the 
computation. Further grid refinement (factor of 3 over the grid used here) may be required to 
achieve a grid independent simulation. The data presented here represent only a small fraction of 
the information contained in the simulation database. 

The results from the direct simulation show no evidence of negative values of turbulent nor- 
mal heat flux in the transitional or turbulent regions. The experimentally observed negative values 
still remain a mystery; however, the present results indicate that it is related to the measurement 
technique. 

In terms of future work, a computation on a refined grid is currently in progress, as are 
detailed flow visualization efforts. A database is being generated that can be used for developing 
transitional models for Reynolds-averaged Navier-Stokes computations, and for developing sub- 
grid scale models for large eddy simulations. 
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